allometry

https://onlinelibrary.wiley.com/doi/abs/10.1111/jeb.12986

explore <- function(y_col, x_col, treatment_col, dt){
  
  log_dt <- data.table(
    y = dt[, get(y_col)],
    x = dt[, get(x_col)],
    treatment = dt[, get(treatment_col)]
  )
  
  fit1 <- lm(y ~ x + treatment,
             data = log_dt,
            na.action = "na.exclude")
  fit2 <- lm(log(y) ~ log(x) + treatment,
             data = log_dt,
            na.action = "na.exclude")
  fit3 <- lm(y/x*100 ~ treatment,
             data = log_dt,
            na.action = "na.exclude")

  gg1 <- ggplot(data = log_dt,
               aes(x = x,
                   y = y,
                   color = treatment)) +
    geom_point() +
    geom_smooth(method = "lm",
                se = FALSE) +
    ylab("muscle weight") +
    xlab("body weight") +
    theme_pubr() +
    theme(legend.title = element_blank()) +
    scale_color_manual(values = pal_okabe_ito_blue)
  
  gg1_log <- ggplot(data = log_dt,
               aes(x = log(x),
                   y = log(y),
                   color = treatment)) +
    geom_point() +
    geom_smooth(method = "lm",
                se = FALSE) +
    ylab("log(muscle weight)") +
    xlab("log(body weight)") +
    theme_pubr() +
    theme(legend.title = element_blank()) +
    scale_color_manual(values = pal_okabe_ito_blue)

  gg2 <- ggplot(data = log_dt,
               aes(x = x,
                   y = y,
                   color = treatment)) +
    geom_point() +
    geom_smooth(method = "lm",
                mapping = aes(y = predict(fit1)),
                se = FALSE) +
    ylab("muscle weight") +
    xlab("body weight") +
    theme_pubr() +
    theme(legend.title = element_blank()) +
    scale_color_manual(values = pal_okabe_ito_blue)

  gg2_log <- ggplot(data = log_dt,
               aes(x = log(x),
                   y = log(y),
                   color = treatment)) +
    geom_point() +
    geom_smooth(method = "lm",
                mapping = aes(y = predict(fit2)),
                se = FALSE) +
    ylab("log(muscle weight)") +
    xlab("log(body weight)") +
    theme_pubr() +
    theme(legend.title = element_blank()) +
    scale_color_manual(values = pal_okabe_ito_blue)

  gg3 <- ggplot(data = log_dt,
               aes(x = x,
                   y = y/x*100,
                   color = treatment)) +
    geom_point() +
    geom_smooth(method = "lm",
                se = FALSE) +
    ylab("Percent Muscle Weight") +
    xlab("body weight") +
    theme_pubr() +
    theme(legend.title = element_blank()) +
    scale_color_manual(values = pal_okabe_ito_blue)
  
  gg <- plot_grid(gg1, gg2, gg3, ncol = 3)
#  gg <- plot_grid(gg1, gg2, gg3, gg4, nrow = 2)
  return(gg)
}

1 ymr

tl;dr – slope 0.149 (-0.137, 0.435)

mice: Male C57BL/6J mice purchased from The Jackson Laboratory (Bar Harbor, ME)

At 60 weeks of age, mice were fasted for 4h, and plasma and tissues collected.

age: 5 months

ymr_folder <- "data - projects/ymr"

fn <- "YMR start020717MMCRIv2.xlsx"
file_path <- here(ymr_folder, fn)

# separately read each genotype
data_range <- "a1:m17"
ymr <- read_excel(file_path,
                    sheet = "Sac020518",
                    range = data_range,
                    col_names = TRUE) %>%
  clean_names() %>%
  data.table() 

ymr[, diet := factor(diet, c("CF", "MR"))]

ymr[, quad := muscle_r_l_quad_g]

all_dt <- data.table(NULL)

treatment_col <- "diet"
groups <- unique(ymr[, get(treatment_col)])
muscles <- c("quad")
body <- "bw_02_05_18"
data_set <- "ymr"
for(group_i in groups){
  data_i <- ymr[diet == group_i] # modify this
  for(muscle_j in muscles){
    form_j <- formula(paste("log(",muscle_j, ") ~ log(", body, ")"))
    m1 <- lm(form_j, data = data_i)
    m1_coef <- cbind(coef(summary(m1)),
                     confint(m1))
    all_dt <- rbind(all_dt,
                    data.table(
                      data = data_set,
                      treatment = ifelse(group_i %in% 
                                           c("CF", "Control","WT", "wt", "wt;wt"),
                                         "Control",
                                         "Treated"),
                      group = group_i,
                      n = summary(m1)$df[2] + 2,
                      muscle = muscle_j,
                      t(m1_coef[2,]
                      ))) 
    
  }
}

1.1 Exploration

gg <- explore(y_col = "muscle_r_l_quad_g",
              x_col = "bw_02_05_18",
              treatment_col = "diet",
              ymr)
gg

1.2 Static allometry

1.3 Additive

m1_scat <- lm(log(scat1_g_leg) ~ log(bw_02_05_18) + diet, data = ymr)
m1_pgat <- lm(log(pgat_g) ~ log(bw_02_05_18) + diet, data = ymr)
m1_muscle <- lm(log(muscle_r_l_quad_g) ~ log(bw_02_05_18) + diet, data = ymr)
m1_liver <- lm(log(liver_g) ~ log(bw_02_05_18) + diet, data = ymr)

scat_coef <- cbind(coef(summary(m1_scat)),
                   confint(m1_scat))
pgat_coef <- cbind(coef(summary(m1_pgat)),
                   confint(m1_pgat))
muscle_coef <- cbind(coef(summary(m1_muscle)),
                     confint(m1_muscle))
liver_coef <- cbind(coef(summary(m1_liver)),
                    confint(m1_liver))

table_dt <- rbind(scat_coef[2,], pgat_coef[2,], muscle_coef[2,], liver_coef[2,])

table_dt %>%
  kable(digits = c(3, 3, 1, 4, 3, 3), caption = "Slopes") %>%
  kable_styling() %>%
  pack_rows("scat", 1, 1) %>%
  pack_rows("pgat", 2, 2) %>%
  pack_rows("muscle", 3, 3) %>%
  pack_rows("liver", 4, 4)

Slopes

Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
scat
3.299	0.379	8.7	0.0000	2.480	4.119
pgat
2.341	0.318	7.4	0.0000	1.655	3.027
muscle
0.149	0.133	1.1	0.2812	-0.137	0.435
liver
1.329	0.165	8.1	0.0000	0.973	1.686

1.4 control diet only

m1_scat <- lm(log(scat1_g_leg) ~ log(bw_02_05_18),
              data = ymr[diet == "CF"])
m1_pgat <- lm(log(pgat_g) ~ log(bw_02_05_18),
              data = ymr[diet == "CF"])
m1_muscle <- lm(log(muscle_r_l_quad_g) ~ log(bw_02_05_18),
              data = ymr[diet == "CF"])
m1_liver <- lm(log(liver_g) ~ log(bw_02_05_18),
              data = ymr[diet == "CF"])

scat_coef <- cbind(coef(summary(m1_scat)),
                   confint(m1_scat))
pgat_coef <- cbind(coef(summary(m1_pgat)),
                   confint(m1_pgat))
muscle_coef <- cbind(coef(summary(m1_muscle)),
                     confint(m1_muscle))
liver_coef <- cbind(coef(summary(m1_liver)),
                    confint(m1_liver))

table_dt <- rbind(scat_coef[2,], pgat_coef[2,], muscle_coef[2,], liver_coef[2,])

table_dt %>%
  kable(digits = c(3, 3, 1, 4, 3, 3), caption = "Slopes") %>%
  kable_styling() %>%
  pack_rows("scat", 1, 1) %>%
  pack_rows("pgat", 2, 2) %>%
  pack_rows("muscle", 3, 3) %>%
  pack_rows("liver", 4, 4)

Slopes

Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
scat
3.487	0.294	11.9	0.0000	2.768	4.206
pgat
2.392	0.299	8.0	0.0002	1.662	3.123
muscle
0.095	0.184	0.5	0.6251	-0.356	0.546
liver
1.456	0.197	7.4	0.0003	0.975	1.938

2 Estimates of error variation in Body Weight

2.1 data_01 – Activation of GCN2/ATF4 signals in amygdalar PKC-δ neurons promotes WAT browning under leucine deprivation

What: Four body weight measures per mouse measured over four days gives estimate of error variance due to fluctations in water levels and gut content. This can be used as an “attenuation correction” for estimate of static allometry (slope of log(organ weight) on log(body weight))

Estimated de_attenuation_index: 1.188089

data_from <- "Activation of GCN2-ATF4 signals in amygdalar PKC-δ neurons promotes WAT browning under leucine deprivation"
file_name <- "41467_2020_16662_MOESM4_ESM.xls"
file_path <- here(data_folder, data_from, file_name)

treatment_vec = rep(c("Control", "Leucine-"), each = 6)

weight <- read_excel(file_path,
                     sheet = "Supplementary Figure 2",
                     range = "K5:N18",
                     col_names = TRUE) %>%
  na.omit() %>%
  data.table()
weight[, id := .I]
weight[, treatment := treatment_vec]
weight <- weight %>%
  melt(id.vars = c("treatment", "id"),
       measure.vars = as.character(0:3),
       variable.name = "day",
       value.name = "weight")

fat <- read_excel(file_path,
                     sheet = "Supplementary Figure 2",
                     range = "Q5:T18",
                     col_names = TRUE) %>%
  na.omit() %>%
  data.table()
fat[, id := .I]
fat[, treatment := treatment_vec]
fat <- fat %>%
  melt(id.vars = c("treatment", "id"),
       measure.vars = as.character(0:3),
       variable.name = "day",
       value.name = "fat")

lean <- read_excel(file_path,
                     sheet = "Supplementary Figure 2",
                     range = "W5:Z18",
                     col_names = TRUE) %>%
  na.omit() %>%
  data.table()
lean[, id := .I]
lean[, treatment := treatment_vec]
lean <- lean %>%
  melt(id.vars = c("treatment", "id"),
       measure.vars = as.character(0:3),
       variable.name = "day",
       value.name = "lean")

data_01 <- cbind(weight, fat = fat[, fat], lean = lean[, lean])
data_01[, day := as.numeric(as.character(day))]
data_01[, residual := weight - fat - lean]

data_01_means <- data_01[, .(weight = mean(weight),
                             fat = mean(fat),
                             lean = mean(lean),
                             weight_sd = sd(weight)),
                         by = .(treatment, id)]
data_01_means[treatment == "Control"] %>%
  kable(digits = c(1,1,1,2,1,3),
        caption = "Means and standard deviation of individual mouse weight measured once per day over four days") %>%
  kable_styling()

Means and standard deviation of individual mouse weight measured once per day over four days

treatment	id	weight	fat	lean	weight_sd
Control	1	29.6	3.24	23.9	0.387
Control	2	26.7	3.46	22.4	0.567
Control	3	28.1	2.90	24.5	0.412
Control	4	30.1	2.66	25.5	0.606
Control	5	27.6	3.16	23.3	0.520
Control	6	27.8	3.39	23.9	0.458

m1 <- lmer(weight ~ 1 + (1 | id), data = data_01[treatment == "Control"])
varcor <- VarCorr(m1) %>%
  data.frame %>%
  data.table

sdcor <- varcor[, sdcor]
sd_table <- data.table(
  level = c("among mice", "within mice"),
  sd = sdcor
)
sd_table %>%
  kable(digits = 3,
        caption = "Modeled standard deviations") %>%
  kable_styling(full_width = FALSE)

Modeled standard deviations

level	sd
among mice	1.252
within mice	0.498

atten_table <- data.table(
  measure = c("body weight variance as percent of total",
    "de-attenuation index"),
  value = c(sdcor[2]^2/(sdcor[2]^2 + sdcor[1]^2),
            sdcor[1]^2/(sdcor[1]^2 - sdcor[2]^2))
)

 atten_table %>%
  kable(digits = 3,
        caption = "De-attenuation Index") %>%
  kable_styling(full_width = FALSE)

De-attenuation Index

measure	value
body weight variance as percent of total	0.137
de-attenuation index	1.188

# carroll & ruppert 96
# https://www.tandfonline.com/doi/abs/10.1080/00031305.1996.10473533
# multiply this by b_ols

3 Muscle scaling

3.1 data03 – Chronic muscle weakness and mitochondrial dysfunction in the absence of sustained atrophy in a preclinical sepsis model

tl;dr –

ta 0.337 (-0.277, 0.951) edl 0.472 (-0.487, 1.430) gast 0.069 (-0.557, 0.695) sol 0.183 (-0.563, 0.930)

data: Non-sepsis control (NSC) male data from Fig 3.

mice: Late-middle-aged adult C57BL/6 mice were acquired from the National Institute on Aging, and all experiments were initiated when animals were 16 months old (average male body weight ~34 grams, female body weight ~28 grams).

age: 16 months

# male data from Fig 3d
# female data from supplement Fig

data_from <- "Chronic muscle weakness and mitochondrial dysfunction in the absence of sustained atrophy in a preclinical sepsis model"
file_name <- "elife-49920-fig3-data1-v1.xlsx"
file_path <- here(data_folder, data_from, file_name)

# males
muscle_g <- read_excel(file_path,
                     sheet = "Fig. 3D",
                     range = "A1:E18",
                     col_names = TRUE) %>%
  data.table() %>%
  clean_names()
setnames(muscle_g, "ga", "gast")

muscle_percent <- read_excel(file_path,
                     sheet = "Fig. 3E",
                     range = "A1:E17",
                     col_names = TRUE) %>%
  data.table() %>%
  clean_names()
setnames(muscle_percent, "ga", "gast")

muscle_list <- c("ta", "edl", "gast", "sol")

#find out missing id in percent data
muscle_g[, ratio := as.character(round(edl/ta, 4))]
muscle_percent[, ratio := as.character(round(edl/ta, 4))]
muscle_g[, ratio2 := as.character(round(gast/sol, 4))]
muscle_percent[, ratio2 := as.character(round(gast/sol, 4))]
muscle_g[, ratio3 := as.character(round(edl/sol, 4))]
muscle_percent[, ratio3 := as.character(round(edl/sol, 4))]
muscle_g[, ratio4 := as.character(round(ta/sol, 4))]
muscle_percent[, ratio4 := as.character(round(ta/sol, 4))]

muscle1 <- merge(muscle_g, muscle_percent, by = "ratio")
muscle2 <- merge(muscle_g, muscle_percent, by = "ratio2")
muscle3 <- merge(muscle_g, muscle_percent, by = "ratio3")
muscle4 <- merge(muscle_g, muscle_percent, by = "ratio4")
# dim(muscle1)
# dim(muscle2)
# dim(muscle3)
# dim(muscle4)
male_g <- muscle2
setnames(male_g,
        old = c("ta.x", "edl.x", "gast.x", "sol.x"),
        new = c("ta", "edl", "gast", "sol"))
setnames(male_g,
        old = c("ta.y", "edl.y", "gast.y", "sol.y"),
        new = c("ta_p", "edl_p", "gast_p", "sol_p"))
male_g[, weight_ta := ta/ta_p]
male_g[, weight_edl := edl/edl_p]
male_g[, weight_ga := gast/gast_p]
male_g[, weight_sol := sol/sol_p]
male_g[, body := ta/ta_p]
ycols <- c("ta", "edl", "gast", "sol", "body")
male_g <- male_g[, .SD, .SDcols = ycols]
male_g[, sex := "male"]

# females
# weights are in mg
female_g <- read_excel(file_path,
                     sheet = "Fig. 3-Supp 1B",
                     range = "A1:E8",
                     col_names = TRUE) %>%
  data.table() %>%
  clean_names()
setnames(female_g, "ga", "gast")

# weights are in mg/g
female_fraction <- read_excel(file_path,
                     sheet = "Fig. 3-Supp 1C",
                     range = "A1:E8",
                     col_names = TRUE) %>%
  data.table() %>%
  clean_names()
setnames(female_fraction, "ga", "gast")

# compute body weights based on fractions
female_g[, body_ga := gast/female_fraction[, gast]]
female_g[, body_ta := ta/female_fraction[, ta]]
female_g[, body_edl := edl/female_fraction[, edl]]
female_g[, body_sol := sol/female_fraction[, sol]]
female_g[, body := sol/female_fraction[, sol]]

female_g <- female_g[, .SD, .SDcols = ycols]
female_g[, sex := "female"]

data_03 <- rbind(male_g, female_g)
# convert muscle weight to g from mg
data_03[, gast := gast/1000]
data_03[, ta := ta/1000]
data_03[, edl := edl/1000]
data_03[, sol := sol/1000]

treatment_col <- "sex"
groups <- unique(data_03[, get(treatment_col)])
muscles <- c("ta", "edl", "gast", "sol")
body <- "body"
data_set <- "data_03"
for(group_i in groups){
  data_i <- data_03[sex == group_i] # modify this
  for(muscle_j in muscles){
    form_j <- formula(paste("log(",muscle_j, ") ~ log(", body, ")"))
    m1 <- lm(form_j, data = data_i)
    m1_coef <- cbind(coef(summary(m1)),
                     confint(m1))
    all_dt <- rbind(all_dt,
                    data.table(
                      data = data_set,
                      treatment = ifelse(group_i %in% 
                                           c("CF", "Control","WT", "wt", "wt;wt"),
                                         "Control",
                                         "Treated"),
                      group = group_i,
                      n = summary(m1)$df[2] + 2,
                      muscle = muscle_j,
                      t(m1_coef[2,]
                      ))) 
    
  }
}

3.1.1 Exploration

tibialis anterior

gg <- explore(y_col = "ta",
              x_col = "body",
              treatment_col = "sex",
              data_03)
gg

extensor digitorum longus

gg <- explore(y_col = "edl",
              x_col = "body",
              treatment_col = "sex",
              data_03)
gg

gastrocnemius

gg <- explore(y_col = "gast",
              x_col = "body",
              treatment_col = "sex",
              data_03)
gg

soleus

gg <- explore(y_col = "sol",
              x_col = "body",
              treatment_col = "sex",
              data_03)
gg

3.1.2 Static allometry

3.1.2.1 Additive

m1_ta <- lm(log(ta) ~ log(body) + sex, data = data_03)
ta_coef <- cbind(coef(summary(m1_ta)),
      confint(m1_ta))

m1_edl <- lm(log(edl) ~ log(body) + sex, data = data_03)
edl_coef <- cbind(coef(summary(m1_edl)),
      confint(m1_edl))

m1_ga <- lm(log(gast) ~ log(body) + sex, data = data_03)
ga_coef <- cbind(coef(summary(m1_ga)),
      confint(m1_ga))

m1_sol <- lm(log(sol) ~ log(body) + sex, data = data_03)
sol_coef <- cbind(coef(summary(m1_sol)),
      confint(m1_sol))

table_dt <- rbind(ta_coef[2,], edl_coef[2,], ga_coef[2,], sol_coef[2,])

table_dt %>%
  kable(digits = c(3, 3, 1, 4, 3, 3), caption = "Slopes") %>%
  kable_styling() %>%
  pack_rows("ta", 1, 1) %>%
  pack_rows("edl", 2, 2) %>%
  pack_rows("gast", 3, 3) %>%
  pack_rows("sol", 4, 4)

Slopes

Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
ta
0.052	0.243	0.2	0.8322	-0.454	0.559
edl
0.890	0.503	1.8	0.0920	-0.159	1.939
gast
0.098	0.207	0.5	0.6414	-0.334	0.531
sol
0.308	0.280	1.1	0.2846	-0.276	0.891

3.1.2.2 by treatment

m1_ta_male <- lm(log(ta) ~ log(body), data = data_03[sex == "male"])
ta_male <- cbind(coef(summary(m1_ta_male)),
      confint(m1_ta_male))
m1_ta_female <- lm(log(ta) ~ log(body), data = data_03[sex == "female"])
ta_female <- cbind(coef(summary(m1_ta_female)),
      confint(m1_ta_female))

m1_edl_male <- lm(log(edl) ~ log(body), data = data_03[sex == "male"])
edl_male <- cbind(coef(summary(m1_edl_male)),
      confint(m1_edl_male))
m1_edl_female <- lm(log(edl) ~ log(body), data = data_03[sex == "female"])
edl_female <- cbind(coef(summary(m1_edl_female)),
      confint(m1_edl_female))

m1_ga_male <- lm(log(gast) ~ log(body), data = data_03[sex == "male"])
ga_male <- cbind(coef(summary(m1_ga_male)),
      confint(m1_ga_male))
m1_ga_female <- lm(log(gast) ~ log(body), data = data_03[sex == "female"])
ga_female <- cbind(coef(summary(m1_ga_female)),
      confint(m1_ga_female))

m1_sol_male <- lm(log(sol) ~ log(body), data = data_03[sex == "male"])
sol_male <- cbind(coef(summary(m1_sol_male)),
      confint(m1_sol_male))
m1_sol_female <- lm(log(sol) ~ log(body), data = data_03[sex == "female"])
sol_female <- cbind(coef(summary(m1_sol_female)),
      confint(m1_sol_female))

table_dt <- rbind(ta_male[2,],
                  edl_male[2,],
                  ga_male[2,],
                  sol_male[2,],
                  ta_female[2,],
                  edl_female[2,],
                  ga_female[2,],
                  sol_female[2,])

table_dt %>%
  kable(digits = c(3, 3, 1, 4, 3, 3), caption = "Slopes") %>%
  kable_styling() %>%
  pack_rows("ta - male", 1, 1) %>%
  pack_rows("edl - male", 2, 2) %>%
  pack_rows("gast - male", 3, 3) %>%
  pack_rows("sol - male", 4, 4) %>%
  pack_rows("ta - female", 5, 5) %>%
  pack_rows("edl - female", 6, 7) %>%
  pack_rows("gast - female", 7, 7) %>%
  pack_rows("sol - female", 8, 8)

Slopes

Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
ta - male
0.337	0.286	1.2	0.2586	-0.277	0.951
edl - male
0.472	0.447	1.1	0.3091	-0.487	1.430
gast - male
0.069	0.292	0.2	0.8166	-0.557	0.695
sol - male
0.183	0.348	0.5	0.6069	-0.563	0.930
ta - female
-0.624	0.348	-1.8	0.1325	-1.518	0.269
edl - female
1.883	1.344	1.4	0.2201	-1.572	5.338
gast - female
0.167	0.115	1.5	0.2061	-0.129	0.463
sol - female
0.602	0.475	1.3	0.2603	-0.618	1.823

3.2 data04 – Antagonistic control of myofiber size and muscle protein quality control by the ubiquitin ligase UBR4 during aging

tl;dr – high coef (> 1) in 6 mo, low coef (.1, .4) in 24 mo

ta 6 mos. 1.197 (0.776, 1.618) sol 6 mos. 1.462 (1.120, 1.804) ta 24 mos. 0.114 (-0.397, 0.624) sol 24 mos. 0.417 (0.043, 0.791)

mice: For knockout, UBR4 floxed mice were bred together with ACTA1-Cre mice23 to yield homozygous UBR4fl/fl alleles either with (UBR4 mKO) or without ACTA1- Cre (wild-type controls). At 3 months of age, mice received daily intraperitoneal injections of 1 mg/kg of tamoxifen for 5 days to induce Cre-mediated recombination (UBR4 mKO mice). The wild-type siblings were used as matched controls that were injected with tamoxifen but where no recombination occurred due to the lack of the Cre recombinase.

age: 6 months, 24 months

# data from Fig 6

data_from <- "Antagonistic control of myofiber size and muscle protein quality control by the ubiquitin ligase UBR4 during aging"
file_name <- "41467_2021_21738_MOESM11_ESM.xlsx"
file_path <- here(data_folder, data_from, file_name)

wt <- read_excel(file_path,
                     sheet = "Figure 1i",
                     range = "B2:L3",
                     col_names = FALSE) %>%
  data.table() %>%
  transpose()
wt[, treatment := "WT"]

ko <- read_excel(file_path,
                     sheet = "Figure 1i",
                     range = "P2:AC3",
                     col_names = FALSE) %>%
  data.table() %>%
  transpose()
ko[, treatment := "KO"]

weight_wide <- rbind(wt,ko)
colnames(weight_wide) <- c("6 m", "24 m", "treatment")
weight <- melt(weight_wide,
               id.vars = "treatment",
               measure.vars = c("6 m", "24 m"),
               variable.name = "age",
               value.name = "weight")

## tibialis anterior muscle
wt <- read_excel(file_path,
                     sheet = "Figure 1j",
                     range = "B2:L3",
                     col_names = FALSE) %>%
  data.table() %>%
  transpose()
wt[, treatment := "WT"]

ko <- read_excel(file_path,
                     sheet = "Figure 1j",
                     range = "M2:Z3",
                     col_names = FALSE) %>%
  data.table() %>%
  transpose()
ko[, treatment := "KO"]

ta_wide <- rbind(wt,ko)
colnames(ta_wide) <- c("6 m", "24 m", "treatment")
ta <- melt(ta_wide,
               id.vars = "treatment",
               measure.vars = c("6 m", "24 m"),
               variable.name = "age",
               value.name = "ta")

## soleus muscle
wt <- read_excel(file_path,
                     sheet = "Figure 1k",
                     range = "B2:L3",
                     col_names = FALSE) %>%
  data.table() %>%
  transpose()
wt[, treatment := "WT"]

ko <- read_excel(file_path,
                     sheet = "Figure 1k",
                     range = "M2:Z3",
                     col_names = FALSE) %>%
  data.table() %>%
  transpose()
ko[, treatment := "KO"]

sol_wide <- rbind(wt,ko)
colnames(sol_wide) <- c("6 m", "24 m", "treatment")
sol <- melt(sol_wide,
               id.vars = "treatment",
               measure.vars = c("6 m", "24 m"),
               variable.name = "age",
               value.name = "sol")

data_04 <- cbind(weight, ta = ta[, ta], sol = sol[, sol])
# View(data_04)

# data from authors for Fig 6
data_from <- "Antagonistic control of myofiber size and muscle protein quality control by the ubiquitin ligase UBR4 during aging"
file_name <- "ubr4mko weights etc.xlsx"
file_path <- here(data_folder, data_from, file_name)

## soleus muscle
data_04 <- read_excel(file_path,
                     range = "A1:F50",
                     col_names = TRUE) %>%
  data.table() %>%
  clean_names()
data_04[, c("age", "genotype") := tstrsplit(age_geno, "-", fixed=TRUE)]

# convert muscle weights from mg to g
data_04[, ta := ta/1000]
data_04[, sol := sol/1000]

treatment_col <- "age_geno"
groups <- unique(data_04[, get(treatment_col)])
muscles <- c("ta", "sol")
body <- "body"
data_set <- "data_04"
for(group_i in groups){
  data_i <- data_04[age_geno == group_i] # modify this
  for(muscle_j in muscles){
    form_j <- formula(paste("log(",muscle_j, ") ~ log(", body, ")"))
    m1 <- lm(form_j, data = data_i)
    m1_coef <- cbind(coef(summary(m1)),
                     confint(m1))
    all_dt <- rbind(all_dt,
                    data.table(
                      data = data_set,
                      treatment = ifelse(group_i %in% 
                                           c("CF", "Control","WT", "wt", "wt;wt"),
                                         "Control",
                                         "Treated"),
                      group = group_i,
                      n = summary(m1)$df[2] + 2,
                      muscle = muscle_j,
                      t(m1_coef[2,]
                      ))) 
    
  }
}

3.2.1 Exploration

tibialis anterior

gg <- explore(y_col = "ta",
              x_col = "body",
              treatment_col = "age_geno",
              data_04)
gg

soleus

gg <- explore(y_col = "sol",
              x_col = "body",
              treatment_col = "age_geno",
              data_04)
gg

3.2.2 Static allometry

3.2.2.1 Additive but by age

m1_ta <- lm(log(ta) ~ log(body) + genotype, data = data_04[age == "6m"])
ta_coef_6 <- cbind(coef(summary(m1_ta)),
      confint(m1_ta))

m1_ta <- lm(log(ta) ~ log(body) + genotype, data = data_04[age == "24m"])
ta_coef_24 <- cbind(coef(summary(m1_ta)),
      confint(m1_ta))

m1_sol <- lm(log(sol) ~ log(body) + genotype, data = data_04[age == "6m"])
sol_coef_6 <- cbind(coef(summary(m1_sol)),
      confint(m1_sol))

m1_sol <- lm(log(sol) ~ log(body) + genotype, data = data_04[age == "24m"])
sol_coef_24 <- cbind(coef(summary(m1_sol)),
      confint(m1_sol))

table_dt <- rbind(ta_coef_6[2,], sol_coef_6[2,], ta_coef_24[2,], sol_coef_24[2,])

table_dt %>%
  kable(digits = c(3, 3, 1, 4, 3, 3), caption = "Slopes") %>%
  kable_styling() %>%
  pack_rows("ta 6 mos.", 1, 1) %>%
  pack_rows("sol 6 mos.", 2, 2) %>%
  pack_rows("ta 24 mos.", 3, 3) %>%
  pack_rows("sol 24 mos.", 4, 4)

Slopes

Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
ta 6 mos.
0.698	0.204	3.4	0.0026	0.273	1.123
sol 6 mos.
1.009	0.150	6.7	0.0000	0.697	1.322
ta 24 mos.
0.244	0.139	1.8	0.0962	-0.048	0.535
sol 24 mos.
0.459	0.180	2.6	0.0200	0.081	0.837

3.2.3 Compare inference between ANCOVA and ratio

m2_ta <- lm(ta/body*100 ~ genotype, data = data_04[age == "6m"])
m2_ta_coef_6 <- cbind(coef(summary(m2_ta)),
      confint(m2_ta))

m2_ta <- lm(ta/body*100 ~ genotype, data = data_04[age == "24m"])
m2_ta_coef_24 <- cbind(coef(summary(m2_ta)),
      confint(m2_ta))

m2_sol <- lm(sol/body*100 ~ genotype, data = data_04[age == "6m"])
m2_sol_coef_6 <- cbind(coef(summary(m2_sol)),
      confint(m2_sol))

m2_sol <- lm(sol/body*100 ~ genotype, data = data_04[age == "24m"])
m2_sol_coef_24 <- cbind(coef(summary(m2_sol)),
      confint(m2_sol))

table_dt <- rbind(ta_coef_6[3,], m2_ta_coef_6[2,], sol_coef_6[3,], m2_sol_coef_6[2,], ta_coef_24[3,], m2_ta_coef_24[2,], sol_coef_24[3,] , m2_sol_coef_24[2,])

table_dt %>%
  kable(digits = c(3, 3, 1, 4, 3, 3), caption = "Slopes") %>%
  kable_styling() %>%
  pack_rows("ta 6 mo - ancova", 1, 1) %>%
  pack_rows("ta 6 mo - ratio", 2, 2) %>%
  pack_rows("sol 6 mo - ancova", 3, 3) %>%
  pack_rows("sol 6 mo - ratio", 4, 4) %>%
  pack_rows("ta 24 mo - ancova", 5, 5) %>%
  pack_rows("ta 24 mo - ratio", 6, 6) %>%
  pack_rows("sol 24 mo - ancova", 7, 7) %>%
  pack_rows("sol 24 mo - ratio", 8, 8)

Slopes

Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
ta 6 mo - ancova
-0.139	0.036	-3.9	0.0009	-0.213	-0.064
ta 6 mo - ratio
-0.020	0.005	-3.7	0.0012	-0.031	-0.009
sol 6 mo - ancova
-0.126	0.026	-4.8	0.0001	-0.181	-0.071
sol 6 mo - ratio
-0.003	0.001	-6.4	0.0000	-0.004	-0.002
ta 24 mo - ancova
-0.199	0.030	-6.6	0.0000	-0.261	-0.136
ta 24 mo - ratio
-0.033	0.007	-4.6	0.0002	-0.047	-0.018
sol 24 mo - ancova
-0.046	0.038	-1.2	0.2378	-0.126	0.033
sol 24 mo - ratio
-0.001	0.001	-1.5	0.1564	-0.003	0.001

3.3 data05 – VPS39-deficiency observed in type 2 diabetes impairs muscle stem cell differentiation via altered autophagy and epigenetics

“We observed no differences in body weight or body composition between WT and Vps39+/− mice (Supplementary Fig. 4a–c)” [compared both absolute and percent]

tl;dr –

slope = 0.501 (0.18, 0.822)

inference: even less different.

mice: Wild-type (WT) and Vps39+/− mice49, on a C57BL/ 6J background, were maintained under standard housing conditions with a 12-h light/dark cycle at the animal facility, Sahlgrenska Academy, University of Gothenburg

Skeletal muscle from Vps39+/− mice contained lower Vps39 levels compared to wild-type (WT) littermates (Fig. 5a), confirming that they would be a suitable model to study the role of VPS39 in muscle dysregulation. We observed no dif- ferences in body weight or body composition between WT and Vps39+/− mice (Supplementary Fig. 4a–c). We then examined glucose homeostasis in Vps39+/− mice. First, an OGTT was used to examine glucose tolerance in vivo (Fig. 5b).

An OGTT was carried out for 4-5-month old male and female mice fasted for 5 h (n = 18 for Vps39+/−, and n = 16 for WT mice).

age: 4-5 months

# data from supplement Fig 4
# body
data_from <- "VPS39-deficiency observed in type 2 diabetes impairs muscle stem cell differentiation via altered autophagy and epigenetics"
file_name <- "41467_2021_22068_MOESM11_ESM.xlsx"
file_path <- here(data_folder, data_from, file_name)

genotype_levels <- c("WT", "Vps39+/-")
males <- read_excel(file_path,
                     sheet = "Supplementary Figure 4a",
                     range = "B4:C11",
                     col_names = TRUE) %>%
  data.table() %>%
  melt(measure.vars = genotype_levels,
       variable.name = "genotype",
       value.name = "body")
males[, sex := "male"]

females <- read_excel(file_path,
                     sheet = "Supplementary Figure 4a",
                     range = "E4:F12",
                     col_names = TRUE) %>%
  data.table() %>%
  melt(measure.vars = genotype_levels,
       variable.name = "genotype",
       value.name = "body")
females[, sex := "female"]

data_05 <- rbind(males, females)

# tibialis anterior
males <- read_excel(file_path,
                     sheet = "Supplementary Figure 4c",
                     range = "B5:C12",
                     col_names = TRUE) %>%
  data.table() %>%
  melt(measure.vars = genotype_levels,
       variable.name = "genotype",
       value.name = "ta")
males[, sex := "male"]

females <- read_excel(file_path,
                     sheet = "Supplementary Figure 4c",
                     range = "E5:F13",
                     col_names = TRUE) %>%
  data.table() %>%
  melt(measure.vars = genotype_levels,
       variable.name = "genotype",
       value.name = "ta")
females[, sex := "female"]

ta <- rbind(males, females)

data_05[, ta := ta[, ta]]

data_05[, treatment := paste(genotype, sex)]

treatment_col <- "treatment"
groups <- unique(data_05[, get(treatment_col)])
muscles <- c("ta")
body <- "body"
data_set <- "data_05"
for(group_i in groups){
  data_i <- data_05[treatment == group_i] # modify this
  for(muscle_j in muscles){
    form_j <- formula(paste("log(",muscle_j, ") ~ log(", body, ")"))
    m1 <- lm(form_j, data = data_i)
    m1_coef <- cbind(coef(summary(m1)),
                     confint(m1))
    all_dt <- rbind(all_dt,
                    data.table(
                      data = data_set,
                      treatment = ifelse(group_i %in% 
                                           c("CF", "Control","WT", "wt", "wt;wt"),
                                         "Control",
                                         "Treated"),
                      group = group_i,
                      n = summary(m1)$df[2] + 2,
                      muscle = muscle_j,
                      t(m1_coef[2,]
                      ))) 
    
  }
}

3.3.1 Exploration

tibialis anterior

gg <- explore(y_col = "ta",
              x_col = "body",
              treatment_col = "treatment",
              data_05)
gg

3.3.2 Static allometry

m1_ta <- lm(log(ta) ~ log(body) + genotype*sex, data = data_05)
ta_coef <- cbind(coef(summary(m1_ta)),
      confint(m1_ta))

table_dt <- t(ta_coef[2,])

table_dt %>%
  kable(digits = c(3, 3, 1, 4, 3, 3), caption = "Slope, conditional on genotype and sex") %>%
  kable_styling() %>%
  pack_rows("ta", 1, 1)

Slope, conditional on genotype and sex

Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
ta
0.501	0.155	3.2	0.0037	0.18	0.822

3.3.3 Compare inference from ANCOVA and ratios

m2_ta <- lm(ta/body*100 ~ genotype*sex, data = data_05)

m1_coef <- cbind(coef(summary(m1_ta)),
      confint(m1_ta))
m2_coef <- cbind(coef(summary(m2_ta)),
      confint(m2_ta))

table_dt <- rbind(m1_coef[3:5,], m2_coef[2:4,])

table_dt %>%
  kable(digits = c(5, 3, 2, 4, 3, 3), caption = "Model comparison") %>%
  kable_styling() %>%
  pack_rows("m1 - ancova", 1, 3) %>%
  pack_rows("m2 - percent", 4, 6)

Model comparison

	Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
m1 - ancova
genotypeVps39+/-	0.00074	0.054	0.01	0.9893	-0.111	0.113
sexmale	0.12363	0.061	2.04	0.0529	-0.002	0.249
genotypeVps39+/-:sexmale	-0.00412	0.083	-0.05	0.9609	-0.176	0.168
m2 - percent
genotypeVps39+/-	0.01670	0.017	0.97	0.3424	-0.019	0.052
sexmale	0.00517	0.018	0.29	0.7743	-0.032	0.042
genotypeVps39+/-:sexmale	-0.03399	0.025	-1.34	0.1920	-0.086	0.018

3.4 data06 – Adult stem cell deficits drive Slc29a3 disorders in mice

tl;dr – huh

mice:

age:

# data from Fig 9
# body
data_from <- "Adult stem cell deficits drive Slc29a3 disorders in mice"
file_name <- "41467_2019_10925_MOESM6_ESM.xlsx"
file_path <- here(data_folder, data_from, file_name)

factor_levels <- c("+/+", "-/-", "-/- + AICAR", "-/- + SCT")
data_06 <- read_excel(file_path,
                     sheet = "Fig.9",
                     range = "C5:F11",
                     col_names = TRUE) %>%
  data.table() %>%
  melt(measure.vars = factor_levels,
       variable.name = "treatment",
       value.name = "body")

muscle <- read_excel(file_path,
                     sheet = "Fig.9",
                     range = "C26:AA27",
                     col_names = FALSE) %>%
  data.table() %>%
  transpose(make.names = 1)

colnames(muscle) <- c("sol", "gast")

data_06[, sol := muscle[, sol]]
data_06[, gast := muscle[, gast]]

treatment_col <- "treatment"
groups <- unique(data_06[, get(treatment_col)])
muscles <- c("sol", "gast")
body <- "body"
data_set <- "data_06"
for(group_i in groups){
  data_i <- data_06[treatment == group_i] # modify this
  for(muscle_j in muscles){
    form_j <- formula(paste("log(",muscle_j, ") ~ log(", body, ")"))
    m1 <- lm(form_j, data = data_i)
    m1_coef <- cbind(coef(summary(m1)),
                     confint(m1))
    all_dt <- rbind(all_dt,
                    data.table(
                      data = data_set,
                      treatment = ifelse(group_i %in% 
                                           c("CF", "Control","WT", "wt", "wt;wt"),
                                         "Control",
                                         "Treated"),
                      group = group_i,
                      n = summary(m1)$df[2] + 2,
                      muscle = muscle_j,
                      t(m1_coef[2,]
                      ))) 
    
  }
}

3.4.1 Exploration

soleus

gg <- explore(y_col = "sol",
              x_col = "body",
              treatment_col = "treatment",
              data_06)
gg

gastrocnemius

gg <- explore(y_col = "gast",
              x_col = "body",
              treatment_col = "treatment",
              data_06)
gg

huh

3.5 data07 – Sestrins are evolutionarily conserved mediators of exercise benefits

tl;dr –

mice:

age:

# data from supplement Fig 6
data_from <- "Sestrins are evolutionarily conserved mediators of exercise benefits"
file_name <- "41467_2019_13442_MOESM3_ESM.xlsx"
file_path <- here(data_folder, data_from, file_name)

# WT
wt_sol <- read_excel(file_path,
                     sheet = "Supplementary Figure 6",
                     range = "A29:C34",
                     col_names = TRUE) %>%
  data.table()
colnames(wt_sol) <- c("treatment", "sol_mg", "sol_frac")
wt_sol[, body_sol := sol_mg/1000/sol_frac]

wt_gtn <- read_excel(file_path,
                     sheet = "Supplementary Figure 6",
                     range = "G29:H34",
                     col_names = TRUE) %>%
  data.table()
colnames(wt_gtn) <- c("gast_mg", "gtn_frac")
wt_gtn[, body_gtn := gast_mg/1000/gtn_frac]

wt <- cbind(wt_sol, wt_gtn)

# TKO

tko_sol <- read_excel(file_path,
                     sheet = "Supplementary Figure 6",
                     range = "A40:C43",
                     col_names = FALSE) %>%
  data.table()
colnames(tko_sol) <- c("treatment", "sol_mg", "sol_frac")
tko_sol[, body_sol := sol_mg/1000/sol_frac]

tko_gtn <- read_excel(file_path,
                     sheet = "Supplementary Figure 6",
                     range = "G40:H43",
                     col_names = FALSE) %>%
  data.table()
colnames(tko_gtn) <- c("gast_mg", "gtn_frac")
tko_gtn[, body_gtn := gast_mg/1000/gtn_frac]

tko <- cbind(tko_sol, tko_gtn)

data_07 <- rbind(wt,tko)
data_07[, body := body_sol]
data_07[, sol := sol_mg/1000]
data_07[, gast := gast_mg/1000]

treatment_col <- "treatment"
groups <- unique(data_07[, get(treatment_col)])
muscles <- c("sol", "gast")
body <- "body"
data_set <- "data_07"
for(group_i in groups){
  data_i <- data_07[treatment == group_i] # modify this
  for(muscle_j in muscles){
    form_j <- formula(paste("log(",muscle_j, ") ~ log(", body, ")"))
    m1 <- lm(form_j, data = data_i)
    m1_coef <- cbind(coef(summary(m1)),
                     confint(m1))
    all_dt <- rbind(all_dt,
                    data.table(
                      data = data_set,
                      treatment = ifelse(group_i %in% 
                                           c("CF", "Control","WT", "wt", "wt;wt"),
                                         "Control",
                                         "Treated"),
                      group = group_i,
                      n = summary(m1)$df[2] + 2,
                      muscle = muscle_j,
                      t(m1_coef[2,]
                      ))) 
    
  }
}

3.5.1 Exploration

soleus

gg <- explore(y_col = "sol",
              x_col = "body",
              treatment_col = "treatment",
              data_07)
gg

gastrocnemius

gg <- explore(y_col = "gast",
              x_col = "body",
              treatment_col = "treatment",
              data_07)
gg

3.5.2 static allometry

m1_sol <- lm(log(sol) ~ log(body) + treatment, data = data_07)
sol_coef <- cbind(coef(summary(m1_sol)),
      confint(m1_sol))
m1_gtn <- lm(log(gast) ~ log(body) + treatment, data = data_07)
gtn_coef <- cbind(coef(summary(m1_gtn)),
      confint(m1_gtn))

table_dt <- rbind(sol_coef[2,],
                  gtn_coef[2,])

table_dt %>%
  kable(digits = c(3, 3, 1, 4, 3, 3), caption = "Slopes") %>%
  kable_styling() %>%
  pack_rows("sol", 1, 1) %>%
  pack_rows("gast", 2, 2)

Slopes

Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
sol
1.472	0.638	2.3	0.0604	-0.089	3.033
gast
0.113	0.311	0.4	0.7284	-0.649	0.875

3.5.3 Compare inference from ANCOVA and ratios

m1_sol <- lm(log(sol) ~ log(body) + treatment, data = data_07)
m2_sol <- lm(sol/body*100 ~ treatment, data = data_07)
m1_gtn <- lm(log(gast) ~ log(body) + treatment, data = data_07)
m2_gtn <- lm(gast/body*100 ~ treatment, data = data_07)

m1_sol_coef <- cbind(coef(summary(m1_sol)),
      confint(m1_sol))
m2_sol_coef <- cbind(coef(summary(m2_sol)),
      confint(m2_sol))
m1_gtn_coef <- cbind(coef(summary(m1_gtn)),
      confint(m1_gtn))
m2_gtn_coef <- cbind(coef(summary(m2_gtn)),
      confint(m2_gtn))

table_dt <- rbind(m1_sol_coef[3,],
                  m2_sol_coef[2,],
                  m1_gtn_coef[3,],
                  m2_gtn_coef[2,])

table_dt %>%
  kable(digits = c(5, 3, 2, 4, 3, 3), caption = "Model comparison") %>%
  kable_styling() %>%
  pack_rows("soleus m1 - ancova", 1, 1) %>%
  pack_rows("soleus m2 - percent", 2, 2) %>%
  pack_rows("gastroc m1 - percent", 3, 3) %>%
  pack_rows("gastroc m2 - percent", 4, 4)

Model comparison

Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
soleus m1 - ancova
0.38498	0.133	2.89	0.0279	0.058	0.711
soleus m2 - percent
0.00951	0.003	3.18	0.0156	0.002	0.017
gastroc m1 - percent
0.19331	0.065	2.97	0.0250	0.034	0.353
gastroc m2 - percent
0.11579	0.030	3.81	0.0067	0.044	0.188

3.6 data08 – The neuromuscular junction is a focal point of mTORC1 signaling in sarcopenia

tl;dr –

mice:

age:

# data from supplement Fig 6
data_from <- "The neuromuscular junction is a focal point of mTORC1 signaling in sarcopenia"
file_name <- "41467_2020_18140_MOESM3_ESM.xlsx"
file_path <- here(data_folder, data_from, file_name)

data_08 <- read_excel(file_path,
                     sheet = "Figure 2A, 2B, S2A-B, S2C-F",
                     range = "A2:BI11",
                     col_names = TRUE) %>%
  data.table() %>%
  transpose(make.names = 1,
            keep.names = "label") %>%
  clean_names()

setnames(data_08, "pla", "plant")

treatment_levels <- c("10mCON", "30mCON", "30mRM")
data_08[, treatment := rep(treatment_levels, c(17,20,23))]

treatment_col <- "treatment"
groups <- unique(data_08[, get(treatment_col)])
muscles <- c("ta", "edl", "tri", "quad", "plant", "sol", "gast")
body <- "body_mass"
data_set <- "data_08"
for(group_i in groups){
  data_i <- data_08[treatment == group_i] # modify this
  for(muscle_j in muscles){
    form_j <- formula(paste("log(",muscle_j, ") ~ log(", body, ")"))
    m1 <- lm(form_j, data = data_i)
    m1_coef <- cbind(coef(summary(m1)),
                     confint(m1))
    all_dt <- rbind(all_dt,
                    data.table(
                      data = data_set,
                      treatment = ifelse(group_i %in% 
                                           c("CF", "Control","WT", "wt", "wt;wt"),
                                         "Control",
                                         "Treated"),
                      group = group_i,
                      n = summary(m1)$df[2] + 2,
                      muscle = muscle_j,
                      t(m1_coef[2,]
                      ))) 
    
  }
}

3.6.1 Exploration

tibialis anterior

gg <- explore(y_col = "ta",
              x_col = "body_mass",
              treatment_col = "treatment",
              data_08)
gg

extensor digitorum longus

gg <- explore(y_col = "edl",
              x_col = "body_mass",
              treatment_col = "treatment",
              data_08)
gg

triceps

gg <- explore(y_col = "tri",
              x_col = "body_mass",
              treatment_col = "treatment",
              data_08)
gg

quadruceps

gg <- explore(y_col = "quad",
              x_col = "body_mass",
              treatment_col = "treatment",
              data_08)
gg

plantaris

gg <- explore(y_col = "plant",
              x_col = "body_mass",
              treatment_col = "treatment",
              data_08)
gg

soleus

gg <- explore(y_col = "sol",
              x_col = "body_mass",
              treatment_col = "treatment",
              data_08)
gg

gastrocnemius

gg <- explore(y_col = "gast",
              x_col = "body_mass",
              treatment_col = "treatment",
              data_08)
gg

3.6.2 static allometry

Additive

m1_ta <- lm(log(ta) ~ log(body_mass) + treatment, data = data_08)
m1_edl <- lm(log(edl) ~ log(body_mass) + treatment, data = data_08)
m1_tri <- lm(log(tri) ~ log(body_mass) + treatment, data = data_08)
m1_quad <- lm(log(quad) ~ log(body_mass) + treatment, data = data_08)
m1_pla <- lm(log(plant) ~ log(body_mass) + treatment, data = data_08)
m1_sol <- lm(log(sol) ~ log(body_mass) + treatment, data = data_08)
m1_gast <- lm(log(gast) ~ log(body_mass) + treatment, data = data_08)

m1_ta_coef <- cbind(coef(summary(m1_ta)),
                    confint(m1_ta))
m1_edl_coef <- cbind(coef(summary(m1_edl)),
                    confint(m1_edl))
m1_tri_coef <- cbind(coef(summary(m1_tri)),
                    confint(m1_tri))
m1_quad_coef <- cbind(coef(summary(m1_quad)),
                    confint(m1_quad))
m1_pla_coef <- cbind(coef(summary(m1_pla)),
                    confint(m1_pla))
m1_sol_coef <- cbind(coef(summary(m1_sol)),
                    confint(m1_sol))
m1_gast_coef <- cbind(coef(summary(m1_gast)),
                    confint(m1_gast))

table_dt <- rbind(m1_ta_coef[2,],
                  m1_edl_coef[2,],
                  m1_tri_coef[2,],
                  m1_quad_coef[2,],
                  m1_pla_coef[2,],
                  m1_sol_coef[2,],
                  m1_gast_coef[2,])

table_dt %>%
  kable(digits = c(5, 3, 2, 4, 3, 3), caption = "Model comparison") %>%
  kable_styling() %>%
  pack_rows("ta", 1, 1) %>%
  pack_rows("edl", 2, 2) %>%
  pack_rows("tri", 3, 3) %>%
  pack_rows("quad", 4, 4) %>%
  pack_rows("plant", 5, 5) %>%
  pack_rows("sol", 6, 6) %>%
  pack_rows("gast", 7, 7)

Model comparison

Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
ta
0.63059	0.103	6.12	0e+00	0.424	0.837
edl
0.68807	0.100	6.88	0e+00	0.488	0.888
tri
0.61822	0.093	6.67	0e+00	0.432	0.804
quad
0.56335	0.088	6.43	0e+00	0.388	0.739
plant
0.69801	0.135	5.16	0e+00	0.427	0.969
sol
0.79029	0.212	3.73	5e-04	0.366	1.215
gast
0.62807	0.091	6.92	0e+00	0.446	0.810

Conditional

m2_ta <- lm(log(ta) ~ log(body_mass) * treatment, data = data_08)
m2_edl <- lm(log(edl) ~ log(body_mass) * treatment, data = data_08)
m2_tri <- lm(log(tri) ~ log(body_mass) * treatment, data = data_08)
m2_quad <- lm(log(quad) ~ log(body_mass) * treatment, data = data_08)
m2_pla <- lm(log(plant) ~ log(body_mass) * treatment, data = data_08)
m2_sol <- lm(log(sol) ~ log(body_mass) * treatment, data = data_08)
m2_gast <- lm(log(gast) ~ log(body_mass) * treatment, data = data_08)

m2_ta_coef <- cbind(coef(summary(m2_ta)),
                    confint(m2_ta))
m2_edl_coef <- cbind(coef(summary(m2_edl)),
                    confint(m2_edl))
m2_tri_coef <- cbind(coef(summary(m2_tri)),
                    confint(m2_tri))
m2_quad_coef <- cbind(coef(summary(m2_quad)),
                    confint(m2_quad))
m2_pla_coef <- cbind(coef(summary(m2_pla)),
                    confint(m2_pla))
m2_sol_coef <- cbind(coef(summary(m2_sol)),
                    confint(m2_sol))
m2_gast_coef <- cbind(coef(summary(m2_gast)),
                    confint(m2_gast))

table_dt <- rbind(m2_ta_coef[c(2,5,6),],
                  m2_edl_coef[c(2,5,6),],
                  m2_tri_coef[c(2,5,6),],
                  m2_quad_coef[c(2,5,6),],
                  m2_pla_coef[c(2,5,6),],
                  m2_sol_coef[c(2,5,6),],
                  m2_gast_coef[c(2,5,6),])

table_dt %>%
  kable(digits = c(5, 3, 2, 4, 3, 3), caption = "Model comparison") %>%
  kable_styling() %>%
  pack_rows("ta", 1, 3) %>%
  pack_rows("edl", 4, 6) %>%
  pack_rows("tri", 7, 9) %>%
  pack_rows("quad", 10, 12) %>%
  pack_rows("plant", 13, 15) %>%
  pack_rows("sol", 16, 18) %>%
  pack_rows("gast", 19, 21)

Model comparison

	Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
ta
log(body_mass)	0.58625	0.207	2.84	0.0064	0.172	1.000
log(body_mass):treatment30mCON	0.00498	0.252	0.02	0.9843	-0.500	0.510
log(body_mass):treatment30mRM	0.19197	0.305	0.63	0.5313	-0.419	0.803
edl
log(body_mass)	0.45872	0.191	2.40	0.0199	0.075	0.842
log(body_mass):treatment30mCON	0.16296	0.233	0.70	0.4873	-0.304	0.630
log(body_mass):treatment30mRM	0.66005	0.282	2.34	0.0230	0.095	1.226
tri
log(body_mass)	0.46997	0.185	2.54	0.0141	0.098	0.842
log(body_mass):treatment30mCON	0.20260	0.226	0.90	0.3736	-0.250	0.655
log(body_mass):treatment30mRM	0.19066	0.273	0.70	0.4885	-0.357	0.739
quad
log(body_mass)	0.47532	0.171	2.78	0.0075	0.133	0.818
log(body_mass):treatment30mCON	0.00398	0.208	0.02	0.9848	-0.414	0.422
log(body_mass):treatment30mRM	0.39546	0.252	1.57	0.1226	-0.110	0.901
plant
log(body_mass)	0.58678	0.271	2.17	0.0347	0.044	1.130
log(body_mass):treatment30mCON	0.08885	0.330	0.27	0.7888	-0.573	0.751
log(body_mass):treatment30mRM	0.29629	0.400	0.74	0.4616	-0.505	1.097
sol
log(body_mass)	0.60368	0.426	1.42	0.1622	-0.250	1.458
log(body_mass):treatment30mCON	0.27536	0.519	0.53	0.5979	-0.765	1.316
log(body_mass):treatment30mRM	0.19067	0.628	0.30	0.7627	-1.069	1.450
gast
log(body_mass)	0.56088	0.177	3.17	0.0025	0.206	0.916
log(body_mass):treatment30mCON	-0.03250	0.216	-0.15	0.8808	-0.465	0.400
log(body_mass):treatment30mRM	0.38807	0.261	1.49	0.1431	-0.136	0.912

3.6.3 Compare inference from ANCOVA and ratios

m3_ta <- lm(ta/body_mass*100 ~ treatment, data = data_08)
m3_edl <- lm(edl/body_mass*100 ~ treatment, data = data_08)
m3_tri <- lm(tri/body_mass*100 ~ treatment, data = data_08)
m3_quad <- lm(quad/body_mass*100 ~ treatment, data = data_08)
m3_pla <- lm(plant/body_mass*100 ~ treatment, data = data_08)
m3_sol <- lm(sol/body_mass*100 ~ treatment, data = data_08)
m3_gast <- lm(gast/body_mass*100 ~ treatment, data = data_08)

m3_ta_coef <- cbind(coef(summary(m3_ta)),
                    confint(m3_ta))
m3_edl_coef <- cbind(coef(summary(m3_edl)),
                    confint(m3_edl))
m3_tri_coef <- cbind(coef(summary(m3_tri)),
                    confint(m3_tri))
m3_quad_coef <- cbind(coef(summary(m3_quad)),
                    confint(m3_quad))
m3_pla_coef <- cbind(coef(summary(m3_pla)),
                    confint(m3_pla))
m3_sol_coef <- cbind(coef(summary(m3_sol)),
                    confint(m3_sol))
m3_gast_coef <- cbind(coef(summary(m3_gast)),
                    confint(m3_gast))

table_dt <- rbind(m1_ta_coef[3:4,],
                  m3_ta_coef[2:3,],
                  m1_edl_coef[3:4,],
                  m3_edl_coef[2:3,],
                  m1_tri_coef[3:4,],
                  m3_tri_coef[2:3,],
                  m1_quad_coef[3:4,],
                  m3_quad_coef[2:3,],
                  m1_pla_coef[3:4,],
                  m3_pla_coef[2:3,],
                  m1_sol_coef[3:4,],
                  m3_sol_coef[2:3,],
                  m1_gast_coef[3:4,],
                  m3_gast_coef[2:3,])

table_dt %>%
  kable(digits = c(5, 3, 2, 5, 3, 3), caption = "Model comparison") %>%
  kable_styling() %>%
  pack_rows("ta - ancova", 1, 2) %>%
  pack_rows("ta - ratio", 3, 4) %>%
  pack_rows("edl - ancova", 5, 6) %>%
  pack_rows("edl - ratio", 7, 8) %>%
  pack_rows("tri - ancova", 9, 10) %>%
  pack_rows("tri - ratio", 11, 12) %>%
  pack_rows("quad - ancova", 13, 14) %>%
  pack_rows("quad - ratio", 15, 16) %>%
  pack_rows("plant - ancova", 17, 18) %>%
  pack_rows("plant - ratio", 19, 20) %>%
  pack_rows("sol - ancova", 21, 22) %>%
  pack_rows("sol - ratio", 23, 24) %>%
  pack_rows("gast - ancova", 25, 26) %>%
  pack_rows("gast - ratio", 27, 28)

Model comparison

	Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
ta - ancova
treatment30mCON	-0.13339	0.023	-5.74	0.00000	-0.180	-0.087
treatment30mRM	-0.02563	0.026	-1.00	0.32248	-0.077	0.026
ta - ratio
treatment30mCON	-21.04546	4.505	-4.67	0.00002	-30.066	-12.025
treatment30mRM	3.29310	4.368	0.75	0.45396	-5.453	12.039
edl - ancova
treatment30mCON	-0.11838	0.023	-5.24	0.00000	-0.164	-0.073
treatment30mRM	-0.05747	0.025	-2.31	0.02484	-0.107	-0.008
edl - ratio
treatment30mCON	-3.83175	0.871	-4.40	0.00005	-5.575	-2.088
treatment30mRM	-0.72891	0.844	-0.86	0.39149	-2.419	0.961
tri - ancova
treatment30mCON	-0.17948	0.021	-8.58	0.00000	-0.221	-0.138
treatment30mRM	-0.04729	0.023	-2.05	0.04532	-0.094	-0.001
tri - ratio
treatment30mCON	-73.01226	10.418	-7.01	0.00000	-93.875	-52.150
treatment30mRM	-0.48171	10.101	-0.05	0.96213	-20.709	19.746
quad - ancova
treatment30mCON	-0.25820	0.020	-13.07	0.00000	-0.298	-0.219
treatment30mRM	-0.19600	0.022	-8.98	0.00000	-0.240	-0.152
quad - ratio
treatment30mCON	-163.29789	15.280	-10.69	0.00000	-193.896	-132.700
treatment30mRM	-101.31607	14.815	-6.84	0.00000	-130.983	-71.649
plant - ancova
treatment30mCON	-0.21971	0.030	-7.20	0.00000	-0.281	-0.159
treatment30mRM	-0.17190	0.034	-5.10	0.00000	-0.239	-0.104
plant - ratio
treatment30mCON	-10.72332	1.510	-7.10	0.00000	-13.747	-7.700
treatment30mRM	-7.10006	1.464	-4.85	0.00001	-10.031	-4.169
sol - ancova
treatment30mCON	-0.26485	0.048	-5.54	0.00000	-0.361	-0.169
treatment30mRM	-0.13904	0.053	-2.63	0.01093	-0.245	-0.033
sol - ratio
treatment30mCON	-7.14038	1.339	-5.33	0.00000	-9.821	-4.460
treatment30mRM	-3.34736	1.298	-2.58	0.01251	-5.946	-0.748
gast - ancova
treatment30mCON	-0.24255	0.020	-11.84	0.00000	-0.284	-0.202
treatment30mRM	-0.24329	0.023	-10.76	0.00000	-0.289	-0.198
gast - ratio
treatment30mCON	-96.22884	9.060	-10.62	0.00000	-114.371	-78.086
treatment30mRM	-84.60001	8.784	-9.63	0.00000	-102.190	-67.010

3.7 data09 – Triggering typical nemaline myopathy with compound heterozygous nebulin mutations reveals myofilament structural changes as pathomechanism

tl;dr –

data – Fig 2a, b. Archived data incomplete. Full data from authors.

mice:

age:

# data from supplement Fig 6
data_from <- "Triggering typical nemaline myopathy with compound heterozygous nebulin mutations reveals myofilament structural changes as pathomechanism"
file_name <- "Compound nebulin model tissue weights Jeff.xlsx"
file_path <- here(data_folder, data_from, file_name)

# Neb106KI Male
neb106KI_m <- read_excel(file_path,
                     sheet = "Neb106KI Male",
                     range = "A2:AB32",
                     col_names = TRUE) %>%
  data.table() %>%
  clean_names()

neb106KI_m <- neb106KI_m[c(1:7,13:30)] # dump empty rows

# Neb106KI Female
neb106KI_f <- read_excel(file_path,
                     sheet = "Neb106KI Female",
                     range = "A2:AB35",
                     col_names = TRUE) %>%
  data.table() %>%
  clean_names()

neb106KI_f <- neb106KI_f[c(1:12, 19:33)] # dump empty rows

# Neb106;Ex55 Male
neb106Ex55_m <- read_excel(file_path,
                     sheet = "Neb106;Ex55 Male",
                     range = "A2:AB36",
                     col_names = TRUE) %>%
  data.table() %>%
  clean_names()

neb106Ex55_m <- neb106Ex55_m[c(1:14,19:34)] # dump empty rows

# Neb106;Ex55 Female
neb106Ex55_f <- read_excel(file_path,
                     sheet = "Neb106;Ex55 Female",
                     range = "A2:AB25",
                     col_names = TRUE) %>%
  data.table() %>%
  clean_names()

neb106Ex55_f <- neb106Ex55_f[c(1:9,14:23)] # dump empty rows

# check if colnames are the same for all four sets
# names(neb106KI_m) == names(neb106KI_f)
# names(neb106KI_m) == names(neb106Ex55_m)
# names(neb106KI_m) == names(neb106Ex55_f)

data_09 <- rbind(neb106KI_m, neb106KI_f, neb106Ex55_m, neb106Ex55_f)

data_09[, genotype := pcr]
# fix highlighted values in "PCR" column (4 total)
data_09[genotype == "het;wt", genotype := "het;het"]
data_09[genotype == "wt;het", genotype := "het;het"]

genotype_levels <- c("wt", "hom", "wt;wt","het;het")
data_09[, genotype := factor(genotype, levels = genotype_levels)]
data_09[, genotype_sex := paste(genotype, sex)]

data_09[, tibia := (l_tibia + r_tibia)/2]
data_09[, tib_cran := (l_tib_cran + r_tib_cran)/2]
data_09[, edl := (l_edl + r_edl)/2]
data_09[, quad := (l_quad + r_quad)/2]
data_09[, sol := (l_sol + r_sol)/2]
data_09[, plant := (l_plant + r_plant)/2]
data_09[, gast := (l_gast + r_gast)/2]

treatment_col <- "genotype_sex"
groups <- unique(data_09[, get(treatment_col)])
muscles <- c("tib_cran", "edl", "quad", "sol", "plant", "gast", "diaph")
body <- "bw"
data_set <- "data_09"
for(group_i in groups){
  data_i <- data_09[genotype_sex == group_i] # modify this
  for(muscle_j in muscles){
    form_j <- formula(paste("log(",muscle_j, ") ~ log(", body, ")"))
    m1 <- lm(form_j, data = data_i)
    m1_coef <- cbind(coef(summary(m1)),
                     confint(m1))
    all_dt <- rbind(all_dt,
                    data.table(
                      data = data_set,
                      treatment = ifelse(group_i %in% 
                                           c("CF", "Control","WT", "wt", "wt;wt"),
                                         "Control",
                                         "Treated"),
                      group = group_i,
                      n = summary(m1)$df[2] + 2,
                      muscle = muscle_j,
                      t(m1_coef[2,]
                      ))) 
    
  }
}

3.7.1 Exploration

tibialis cranalis

gg1 <- explore(y_col = "tib_cran",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "M"])

gg2 <- explore(y_col = "tib_cran",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "F"])
plot_grid(gg1,gg2,nrow=2, labels = c("Males", "Females"))

extensor digitorum longus

gg1 <- explore(y_col = "edl",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "M"])

gg2 <- explore(y_col = "edl",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "F"])
plot_grid(gg1,gg2,nrow=2, labels = c("Males", "Females"))

quadriceps

gg1 <- explore(y_col = "quad",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "M"])

gg2 <- explore(y_col = "quad",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "F"])
plot_grid(gg1,gg2,nrow=2, labels = c("Males", "Females"))

soleus

gg1 <- explore(y_col = "sol",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "M"])

gg2 <- explore(y_col = "sol",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "F"])
plot_grid(gg1,gg2,nrow=2, labels = c("Males", "Females"))

plantaris

gg1 <- explore(y_col = "plant",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "M"])

gg2 <- explore(y_col = "plant",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "F"])
plot_grid(gg1,gg2,nrow=2, labels = c("Males", "Females"))

gastrocnemius

gg1 <- explore(y_col = "gast",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "M"])

gg2 <- explore(y_col = "gast",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "F"])
plot_grid(gg1,gg2,nrow=2, labels = c("Males", "Females"))

diaphragm

gg1 <- explore(y_col = "diaph",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "M"])

gg2 <- explore(y_col = "diaph",
              x_col = "bw",
              treatment_col = "genotype",
              data_09[sex == "F"])
plot_grid(gg1,gg2,nrow=2, labels = c("Males", "Females"))

3.7.2 Static allometry

Additive

m1_tib_cran <- lm(log(tib_cran) ~ log(bw) + genotype_sex, data = data_09)
m1_edl <- lm(log(edl) ~ log(bw) + genotype_sex, data = data_09)
m1_quad <- lm(log(quad) ~ log(bw) + genotype_sex, data = data_09)
m1_pla <- lm(log(plant) ~ log(bw) + genotype_sex, data = data_09)
m1_sol <- lm(log(sol) ~ log(bw) + genotype_sex, data = data_09)
m1_gast <- lm(log(gast) ~ log(bw) + genotype_sex, data = data_09)
m1_diaph <- lm(log(diaph) ~ log(bw) + genotype_sex, data = data_09)

m1_tib_cran_coef <- cbind(coef(summary(m1_tib_cran)),
                    confint(m1_tib_cran))
m1_edl_coef <- cbind(coef(summary(m1_edl)),
                    confint(m1_edl))
m1_quad_coef <- cbind(coef(summary(m1_quad)),
                    confint(m1_quad))
m1_pla_coef <- cbind(coef(summary(m1_pla)),
                    confint(m1_pla))
m1_sol_coef <- cbind(coef(summary(m1_sol)),
                    confint(m1_sol))
m1_gast_coef <- cbind(coef(summary(m1_gast)),
                    confint(m1_gast))
m1_diaph_coef <- cbind(coef(summary(m1_diaph)),
                    confint(m1_diaph))

table_dt <- rbind(m1_tib_cran_coef[2,],
                  m1_edl_coef[2,],
                  m1_quad_coef[2,],
                  m1_pla_coef[2,],
                  m1_sol_coef[2,],
                  m1_gast_coef[2,],
                  m1_diaph_coef[2,]
                  )

table_dt %>%
  kable(digits = c(5, 3, 2, 4, 3, 3), caption = "Model comparison") %>%
  kable_styling() %>%
  pack_rows("tib_cran", 1, 1) %>%
  pack_rows("edl", 2, 2) %>%
  pack_rows("quad", 3, 3) %>%
  pack_rows("plant", 4, 4) %>%
  pack_rows("sol", 5, 5) %>%
  pack_rows("gast", 6, 6) %>%
  pack_rows("diaph", 7, 7)

Model comparison

Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
tib_cran
0.32147	0.072	4.44	0.0000	0.178	0.465
edl
0.29284	0.092	3.19	0.0021	0.110	0.476
quad
0.34569	0.073	4.76	0.0000	0.201	0.490
plant
0.38450	0.080	4.82	0.0000	0.226	0.543
sol
0.53831	0.093	5.81	0.0000	0.354	0.723
gast
0.38785	0.065	6.01	0.0000	0.260	0.516
diaph
0.56971	0.085	6.70	0.0000	0.400	0.739

3.7.3 Inference

males using body weight

m1_tib_cran <- lm(log(tib_cran) ~ log(bw) + genotype,
                  data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m1_edl <- lm(log(edl) ~ log(bw) + genotype,
             data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m1_quad <- lm(log(quad) ~ log(bw) + genotype,
              data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m1_pla <- lm(log(plant) ~ log(bw) + genotype,
             data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m1_sol <- lm(log(sol) ~ log(bw) + genotype,
             data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m1_gast <- lm(log(gast) ~ log(bw) + genotype,
              data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m1_diaph <- lm(log(diaph) ~ log(bw) + genotype,
               data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])

m2_tib_cran <- lm(tib_cran/bw*100 ~ genotype,
                  data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m2_edl <- lm(edl/bw*100 ~ genotype,
             data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m2_quad <- lm(quad/bw*100 ~ genotype,
              data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m2_pla <- lm(plant/bw*100 ~ genotype,
             data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m2_sol <- lm(sol/bw*100 ~ genotype,
             data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m2_gast <- lm(gast/bw*100 ~ genotype,
              data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m2_diaph <- lm(diaph/bw*100 ~ genotype,
               data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])


m1_tib_cran_coef <- cbind(coef(summary(m1_tib_cran)),
                    confint(m1_tib_cran))
m1_edl_coef <- cbind(coef(summary(m1_edl)),
                    confint(m1_edl))
m1_quad_coef <- cbind(coef(summary(m1_quad)),
                    confint(m1_quad))
m1_pla_coef <- cbind(coef(summary(m1_pla)),
                    confint(m1_pla))
m1_sol_coef <- cbind(coef(summary(m1_sol)),
                    confint(m1_sol))
m1_gast_coef <- cbind(coef(summary(m1_gast)),
                    confint(m1_gast))
m1_diaph_coef <- cbind(coef(summary(m1_diaph)),
                    confint(m1_diaph))

m2_tib_cran_coef <- cbind(coef(summary(m2_tib_cran)),
                    confint(m2_tib_cran))
m2_edl_coef <- cbind(coef(summary(m2_edl)),
                    confint(m2_edl))
m2_quad_coef <- cbind(coef(summary(m2_quad)),
                    confint(m2_quad))
m2_pla_coef <- cbind(coef(summary(m2_pla)),
                    confint(m2_pla))
m2_sol_coef <- cbind(coef(summary(m2_sol)),
                    confint(m2_sol))
m2_gast_coef <- cbind(coef(summary(m2_gast)),
                    confint(m2_gast))
m2_diaph_coef <- cbind(coef(summary(m2_diaph)),
                    confint(m2_diaph))

table_dt <- rbind(m1_tib_cran_coef[3,],
                  m2_tib_cran_coef[2,],
                  m1_edl_coef[3,],
                  m2_edl_coef[2,],
                  m1_quad_coef[3,],
                  m2_quad_coef[2,],
                  m1_pla_coef[3,],
                  m2_pla_coef[2,],
                  m1_sol_coef[3,],
                  m2_sol_coef[2,],
                  m1_gast_coef[3,],
                  m2_gast_coef[2,],
                  m1_diaph_coef[3,],
                  m2_diaph_coef[2,]
                  )

table_dt %>%
  kable(digits = c(5, 3, 2, 4, 3, 3), caption = "Model comparison") %>%
  kable_styling() %>%
  pack_rows("tib_cran - ancova", 1, 1) %>%
  pack_rows("tib_cran - ratio", 2, 2) %>%
  pack_rows("edl - ancova", 3, 3) %>%
  pack_rows("edl - ratio", 4, 4) %>%
  pack_rows("quad - ancova", 5, 5) %>%
  pack_rows("quad - ratio", 6, 6) %>%
  pack_rows("plant - ancova", 7, 7) %>%
  pack_rows("plant - ratio", 8, 8) %>%
  pack_rows("sol - ancova", 9, 9) %>%
  pack_rows("sol - ratio", 10, 10) %>%
  pack_rows("gast - ancova", 11, 11) %>%
  pack_rows("gast - ratio", 12, 12) %>%
  pack_rows("diaph - ancova", 13, 13) %>%
  pack_rows("diaph - ratio", 14, 14)

Model comparison

Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
tib_cran - ancova
0.25802	0.039	6.56	0.0000	0.177	0.339
tib_cran - ratio
0.07367	0.007	10.10	0.0000	0.059	0.089
edl - ancova
-0.04227	0.060	-0.70	0.4891	-0.167	0.083
edl - ratio
0.00295	0.002	1.89	0.0722	0.000	0.006
quad - ancova
-0.08423	0.031	-2.73	0.0120	-0.148	-0.020
quad - ratio
0.02490	0.026	0.96	0.3488	-0.029	0.079
plant - ancova
-0.05769	0.044	-1.32	0.2004	-0.148	0.033
plant - ratio
0.00528	0.003	1.94	0.0632	0.000	0.011
sol - ancova
0.10536	0.051	2.06	0.0508	0.000	0.211
sol - ratio
0.00812	0.001	5.60	0.0000	0.005	0.011
gast - ancova
-0.19742	0.030	-6.51	0.0000	-0.260	-0.135
gast - ratio
-0.02419	0.017	-1.42	0.1680	-0.059	0.011
diaph - ancova
0.22855	0.046	4.97	0.0001	0.133	0.324
diaph - ratio
0.11799	0.016	7.39	0.0000	0.085	0.151

males using tibia length (following authors)

m1_tib_cran <- lm(log(tib_cran) ~ log(tibia) + genotype,
                  data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m1_edl <- lm(log(edl) ~ log(tibia) + genotype,
             data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m1_quad <- lm(log(quad) ~ log(tibia) + genotype,
              data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m1_pla <- lm(log(plant) ~ log(tibia) + genotype,
             data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m1_sol <- lm(log(sol) ~ log(tibia) + genotype,
             data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m1_gast <- lm(log(gast) ~ log(tibia) + genotype,
              data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m1_diaph <- lm(log(diaph) ~ log(tibia) + genotype,
               data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])

m2_tib_cran <- lm(tib_cran/tibia*100 ~ genotype,
                  data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m2_edl <- lm(edl/tibia*100 ~ genotype,
             data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m2_quad <- lm(quad/tibia*100 ~ genotype,
              data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m2_pla <- lm(plant/tibia*100 ~ genotype,
             data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m2_sol <- lm(sol/tibia*100 ~ genotype,
             data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m2_gast <- lm(gast/tibia*100 ~ genotype,
              data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])
m2_diaph <- lm(diaph/tibia*100 ~ genotype,
               data = data_09[sex == "M" & genotype %in% c("wt;wt", "het;het")])


m1_tib_cran_coef <- cbind(coef(summary(m1_tib_cran)),
                    confint(m1_tib_cran))
m1_edl_coef <- cbind(coef(summary(m1_edl)),
                    confint(m1_edl))
m1_quad_coef <- cbind(coef(summary(m1_quad)),
                    confint(m1_quad))
m1_pla_coef <- cbind(coef(summary(m1_pla)),
                    confint(m1_pla))
m1_sol_coef <- cbind(coef(summary(m1_sol)),
                    confint(m1_sol))
m1_gast_coef <- cbind(coef(summary(m1_gast)),
                    confint(m1_gast))
m1_diaph_coef <- cbind(coef(summary(m1_diaph)),
                    confint(m1_diaph))

m2_tib_cran_coef <- cbind(coef(summary(m2_tib_cran)),
                    confint(m2_tib_cran))
m2_edl_coef <- cbind(coef(summary(m2_edl)),
                    confint(m2_edl))
m2_quad_coef <- cbind(coef(summary(m2_quad)),
                    confint(m2_quad))
m2_pla_coef <- cbind(coef(summary(m2_pla)),
                    confint(m2_pla))
m2_sol_coef <- cbind(coef(summary(m2_sol)),
                    confint(m2_sol))
m2_gast_coef <- cbind(coef(summary(m2_gast)),
                    confint(m2_gast))
m2_diaph_coef <- cbind(coef(summary(m2_diaph)),
                    confint(m2_diaph))

table_dt <- rbind(m1_tib_cran_coef[3,],
                  m2_tib_cran_coef[2,],
                  m1_edl_coef[3,],
                  m2_edl_coef[2,],
                  m1_quad_coef[3,],
                  m2_quad_coef[2,],
                  m1_pla_coef[3,],
                  m2_pla_coef[2,],
                  m1_sol_coef[3,],
                  m2_sol_coef[2,],
                  m1_gast_coef[3,],
                  m2_gast_coef[2,],
                  m1_diaph_coef[3,],
                  m2_diaph_coef[2,]
                  )

table_dt %>%
  kable(digits = c(5, 3, 2, 4, 3, 3), caption = "Model comparison") %>%
  kable_styling() %>%
  pack_rows("tib_cran - ancova", 1, 1) %>%
  pack_rows("tib_cran - ratio", 2, 2) %>%
  pack_rows("edl - ancova", 3, 3) %>%
  pack_rows("edl - ratio", 4, 4) %>%
  pack_rows("quad - ancova", 5, 5) %>%
  pack_rows("quad - ratio", 6, 6) %>%
  pack_rows("plant - ancova", 7, 7) %>%
  pack_rows("plant - ratio", 8, 8) %>%
  pack_rows("sol - ancova", 9, 9) %>%
  pack_rows("sol - ratio", 10, 10) %>%
  pack_rows("gast - ancova", 11, 11) %>%
  pack_rows("gast - ratio", 12, 12) %>%
  pack_rows("diaph - ancova", 13, 13) %>%
  pack_rows("diaph - ratio", 14, 14)

Model comparison

Estimate	Std. Error	t value	Pr(>|t|)	2.5 %	97.5 %
tib_cran - ancova
0.22766	0.034	6.63	0.0000	0.157	0.298
tib_cran - ratio
0.06878	0.009	7.55	0.0000	0.050	0.088
edl - ancova
-0.12575	0.046	-2.75	0.0118	-0.221	-0.031
edl - ratio
-0.00744	0.002	-3.49	0.0020	-0.012	-0.003
quad - ancova
-0.08897	0.021	-4.30	0.0003	-0.132	-0.046
quad - ratio
-0.11278	0.021	-5.32	0.0000	-0.156	-0.069
plant - ancova
-0.09259	0.038	-2.46	0.0213	-0.170	-0.015
plant - ratio
-0.00889	0.003	-2.74	0.0112	-0.016	-0.002
sol - ancova
0.02933	0.045	0.66	0.5171	-0.063	0.121
sol - ratio
0.00084	0.002	0.40	0.6921	-0.003	0.005
gast - ancova
-0.21667	0.026	-8.22	0.0000	-0.271	-0.162
gast - ratio
-0.15730	0.017	-9.47	0.0000	-0.191	-0.123
diaph - ancova
0.14578	0.047	3.12	0.0048	0.049	0.242
diaph - ratio
0.08380	0.027	3.14	0.0044	0.029	0.139

4 What is the scaling coefficient?

all_dt <- clean_names(all_dt)

4.1 All skeletal muscles combined

4.1.1 Modeled mean slopes

m1 <- lm(estimate ~ 1, data = all_dt)
m1_coef <- cbind(coef(summary(m1)), confint(m1))
m2 <- lm(estimate ~ 1, weights = 1/(all_dt$std_error), data = all_dt)
m2_coef <- cbind(coef(summary(m2)), confint(m2))

set.seed(1)
boot_n <- 2000
mean_i <- numeric(boot_n)
weighted_i <- numeric(boot_n)
full_set <- 1:nrow(all_dt)
inc <- full_set
for(i in 1:boot_n){
  mean_i[i] <- mean(all_dt[inc, estimate])
  weighted_i[i] <- weighted.mean(all_dt[inc, estimate], 1/all_dt[inc, std_error])
  inc <- sample(full_set, replace = TRUE)
}
ci_mean <- quantile(mean_i, c(0.025, 0.975))
ci_weighted <- quantile(weighted_i, c(0.025, 0.975))
out_table <- data.table(
  Method = c("lm", "weighted lm", "boot", "weighted boot"),
  Mean = c(m1_coef[1, "Estimate"], m2_coef[1, "Estimate"], mean_i[1], weighted_i[1]),
  "2.5%" = c(m1_coef[1, "2.5 %"], m2_coef[1, "2.5 %"], ci_mean[1], ci_weighted[1]),
  "97.5%" = c(m1_coef[1, "97.5 %"], m2_coef[1, "97.5 %"], ci_mean[2], ci_weighted[2])
  )

out_table %>%
  kable(digits = c(1,3,3,3)) %>%
  kable_styling()

Method	Mean	2.5%	97.5%
lm	0.439	0.334	0.544
weighted lm	0.452	0.386	0.518
boot	0.439	0.330	0.538
weighted boot	0.452	0.395	0.504

4.2 Individual skeletal muscles

m1 <- lm(estimate ~ 1 + muscle, data = all_dt)
m2 <- lm(estimate ~ 1 + muscle, weights = 1/(all_dt$std_error), data = all_dt)
m3 <- lmer(estimate ~ 1 + (1 | muscle), data = all_dt)

4.2.1 Modeled mean slope for each muscle (weighted means)

emmeans(m2, specs = "muscle") %>%
  kable(digits = c(1,3,3,1,3,3)) %>%
  kable_styling()

muscle	emmean	SE	df	lower.CL	upper.CL
diaph	0.580	0.125	102	0.333	0.828
edl	0.468	0.100	102	0.270	0.666
gast	0.384	0.081	102	0.223	0.544
plant	0.477	0.107	102	0.264	0.689
quad	0.386	0.085	102	0.217	0.555
sol	0.530	0.091	102	0.350	0.711
ta	0.463	0.098	102	0.270	0.657
tib_cran	0.332	0.117	102	0.099	0.564
tri	0.599	0.168	102	0.267	0.931

4.2.2 ANOVA

comb_dt <- data.table(
  " " = c("muscle", "muscle"),
  rbind(anova(m1)[1,], anova(m2)[1,])
)
comb_dt %>%
  kable(digits = c(1, 1,2,2,2,2)) %>%
  kable_styling() %>%
  pack_rows("unweighted", 1, 1) %>%
  pack_rows("weighted", 2, 2)

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
unweighted
muscle	8	0.83	0.10	0.32	0.96
weighted
muscle	8	2.76	0.34	0.63	0.75

4.3 Plot of slope vs. SE

qplot(x = std_error, y = estimate, data = all_dt) +
  geom_hline(yintercept = out_table[Method == "weighted lm", Mean],
             linetype = "dashed") +
  theme_pubr()

m1 <- lm(estimate ~ std_error + muscle, data = all_dt)
qplot(x = std_error, y = estimate, color = muscle, data = all_dt) +
  geom_hline(yintercept = out_table[Method == "weighted lm", Mean],
             linetype = "dashed") +
  theme_pubr()